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ABSTRACT 

The current work illustrates the extension of computational fluid dynamics techniques to artificial heart flow 
simulation. Unsteady incompressible Navier-Stokes equations written in three-dimensional generalized curvilinear 
coordinates are solved iteratively at each physical time step until the incompressibility condition is satisfied. The 
solution method is based on the pseudo-compressibility approach and uses an implicit-upwind differencing 
scheme together with the Gauss-Seidel line relaxation method. The efficiency and robustness of the time-accurate 
formulation of the numerical algorithm are tested by computing the flow through model geometries. A channel 
flow with a moving indentation is computed and validated with experimental measurements and other numerical 
solutions. In order to handle the geometric complexity and the moving boundary problems, a zonal method and 
an overlapped grid embedding scheme are employed, respectively. Steady-state solutions for the flow through 
a tilting-disk heart valve are compared against experimental measurements. Good agreement is obtained. The 
flow computation during the valve opening and closing is carried out to illustrate the moving boundary capability. 
Aided by experimental evidence, the flow through an entire Penn State artificial heart model is computed. 

L INTRODUCTION 

With the advent of supercomputer hardware as well as fast numerical methods, researchers in the field of 
computational fluid dynamics (CFD) tackle more complicated problems than ever before. With these new 
capabilities, CFD has become an essential part of aerospace research and design. For example, the 
incompressible flow solver developed by Kwak et al [1] was extensively used for simulating the flow through space 
shuttle main engine power head components. The redesign of the space shuttle main engine hot gas manifold, 
guided by the computations of Chang et al. [2] illustrates the usefulness of CFD in the aerospace research. 

Extending the CFD technology developed for the aerospace industry to artificial heart simulations will open a 
new route in the artificial heart research. Artificial heart devices have been used widely since the early 1960s to 
replace or to assist natural organs. The replacement can be a heart valve or a total artificial heart. However, the 
prosthetic devices are found to be less efficient than natural organs and various problems have been found during 
clinical observations. The most serious problems are believed to be directly associated with the flow fields created 
inside the artificial heart devices. Major difficulties originating from fluid dynamics phenomena include : 1)- 
Separated and secondary flow regions cause clotting; 2)- High turbulent shear stress can damage the red blood 
cells; 3^Large pressure losses across the valves prevent the heart from working efficiently. Several experimental 
studies on commonly used valve geometries have pointed out the adverse effects of the stagnation and 
recirculation regions on blood flow. Although the experimental studies played an important role in the design 
process of these devices, they can provide flow characteristics for only limited regions of the flow field. In 
addition, the experimental measurements are very difficult because of the moving boundaries in the artificial 
heart devices. Having detailed knowledge of the flow quantities can help a design engineer improve the artificial 
heart and valve geometries, where a smooth flow is desired. The development of such a numerical simulation 
tool, which can be used in the artificial heart research programs, is initiated in the present study. 
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Computational studies of blood flow in hearts and heart valves have been quite limited. The most notable work 
has been performed by Peskin and his co-workers. 6 " 8 Their main effort has been directed to simulations of the 
natural heart and heart valves combining Eulerian flow equations and a Lagrangian description of heart walls 
and valves. Peskin and McQueen 6 modeled the prosthetic heart valves in the numerical simulation of the flow 
in the natural heart. They used boundary forces derived from the energy function in order to model valve 
opening and closing, and they also modeled the elastic behavior of the walls. Their solutions were obtained for 
low Reynolds numbers in two dimensions using a square cartesian mesh. McCracken and Peskin applied a 
combined vortex-grid method for the blood flow through the mitral valve in two dimensions. Peskin and 
McQueen 8 demonstrated the capability of modeling the elastic behavior of heart muscles by applying their 
extended three dimensional solution procedure to a toroidal tube. Although determining the elastic behavior of 
the walls is an important task in the natural heart study, the boundary motion in many artificial heart devices, 
such as the Penn State electric artificial heart, can be well defined. 9 Therefore, the present study is focused on 
the fluid problem with prescribed body motion. 

Underwood and Mueller 10 obtained the flow characteristics for the Kay-Shiley disk type valve using the stream 
function-vortidty formulation. Their results showed agreement with experimental data up to a Reynolds number 
of 600. Idelsohn, Costa, and Ponso 11 modeled the flow through the Kay-Shiley caged disk, Starr-Edwars caged 
ball, and Bjork-Shiley tilting disk valves and compared their performance. Turbulent flow through trileaflet aortic 
heart valves was simulated by Stevensen et al [12]. Most numerical studies assumed that the flow through the 
heart valve was two-dimensional. Additionally, the valve opening and closing motion was neglected; only the flow 
through a fixed valve position was studied. In reality, the geometry is three-dimensional, and the flow through 
heart valves involves moving boundaries. 

The current study proposes the development of a computational procedure simulating steady and unsteady, 
three-dimensional flows through artificial hearts and heart valves with moving boundaries. In the next sections, 
the method of solution is summarized followed by demonstration how this CFD procedure can be used for 
probing the flow through artificial heart devices. 

II, METHOD OF SOLUTION 

The flow through artificial heart devices is unsteady, viscous, and incompressible. In the present study, the 
non-Newtonian nature of the blood is neglected and the flow is described by the three-dimensional 
incompressible Navier-Stokes equations. Numerical simulation of such flows is a very challenging problem in 
computational fluid dynamics. In addition to the geometric complexity, obtaining time-dependent solutions of 
the incompressible Navier-Stokes equations with moving boundaries poses many difficult numerical problems. 

Because the pressure and velocity fields are not directly coupled due to the lack of a pressure term in the 
continuity equation, numerical solution of the incompressible Navier-Stokes equations requires special attention 
in order to satisfy the divergence-free constraint on the velocity field. The most widely used methods which use 
primative variables are fractional-step and pseudocompressibility techniques. In fractional step method, the 
auxiliary velocity field is solved by using the momentum equations. Then, a Poisson equation for pressure is 
formed by taking the divergence of the momentum equations and by using a divergence- free velocity field 
constraint. Solving the Poisson equation for pressure efficiently in three-dimensional curvilinear coordinates is 
the most important feature of the fractional step method. 13 One way to avoid the numerical difficulty originated 
by the elliptic nature of the problem is to use a pseudocompressibility method. With the pseudocompressibility 
method, the elliptic-parabolic type equations are transformed into hyperbolic-parabolic type equations. Well 
established solution algorithms developed for compressible flows can be utilized to solve the resulting equations. 

In the present study, the pseudocompressibility approach is used and the time accuracy is attained by iterating 
in pseudo-time until the divergence of velocity is driven toward zero to within a specified error tolerance. Here, 
the time derivatives in the momentum equations are differenced using a second-order, three-point, backward- 
difference formula. The numerical method uses a second-order central difference for viscous terms and a higher 
order flux-difference splitting for the convective terms. The resulting matrix equation is solved iteratively by using 
a nonfactored line relaxation scheme, which maintains stability and allows a large pseudo-time step to be taken. 
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At each sweep direction, a tridiagonal matrix is formed and off line terms of the matrix equation are moved to 
the right-hand side of the equation. Details of the governing equations and numerical method are given in Refs. 
14-15. 

One of the biggest difficulties in the simulation of flows in complicated three- dimensional configurations is the 
discretization of the physical domain. The problem becomes more severe if one body in the domain of interest 
moves relative to another one as is seen in the tilting disk valve and the Penn State artificial heart geometries. 
The use of a zonal approach is a practical solution if the grids are stationary. For more general applications, a 
Chimera grid embedding technique 16 provides a greater flexibility for the grid motion. This technique is employed 
for flow computations through a tilting disk valve and the Penn State Artificial heart model. 

III. COMPUTED RESULTS 

One of the goals of this study is to simulate the flow through a realistic model of an artificial heart. Since the 
geometry and the flow physics are complicated, the computational procedure is validated by solving several 
simpler problems which characterize the flow in various parts of an artificial heart. As a first step, an idealized 
2-D pump model was chosen to demonstrate the capability of the time-accurate formulation under a moving grid 
condition. Geometry of this model and computed results are presented in Ref. 14. Channel flow with an indented 
wall, the flow through a tilting disk heart valve and the flow through the Penn State artificial heart model are 
included in this section. 


Channel Flow with a Moving Indentation 

Channel flow with an asymmetric oscillating indentation was experimentally studied by Pedley and Stephanoff, 17 
and was numerically simulated by Ralph and Pedley. 18 In the experiment, the channel wall was rigid everywhere 
except for the indentation which is made of a thick rubber membrane. The experiment shows that flow to be 
two-dimensional near the midplane and so the computation is done in two dimensions. 

Figure 1 illustrates the instantaneous streamlines plotted at several nondimensional times for Re = 600 and St 
= 0.057. The Reynolds number is based on the channel height, a, and the average velocity at the entrance of the 
channel, U. Strouhal number is defined as St = af/U where/ is the oscillation frequency. At the beginning of 
the cycle the flow downstream of the indentation is parallel to the channel walls. A single eddy is formed at the 
sloping wall of the indentation during the first half of the cycle. The streamlines at the core flow are lifted slightly 
upwards as shown in Fig. 1-a. This is the beginning of the wavy flow patterns of the core flow. The formation of 
a second separated eddy on the opposite wall can be seell in Fig. I-b. In later stages, the double row of eddies 
along the lower and upper wall of the channel is observed. At the end of the cycle the vortices are swept 
downstream (Fig. 1-g), and the residual vortices are not strong enough to affect the next cycle. In fact, the flow 
patterns of the first and second cycles are quite similar. Consequently, the flow is assumed periodic in time, even 
at the first cycle. Another interesting phenomenon observed in the present study as well as observed in the 
experimental and other numerical studies is the eddy- doubling which can be seen in Figs. 1-c through 1-e. It 
occurred in the second, third, and fourth vortices from the indentation. In eddy-doubling, a single eddy splits into 
two co-rotating eddies. 

The time evolution of the vortices' centers is compared against the experimental and other numerical findings. 
The first four vortices are called Vortices A, B, C, and D, and are shown in Fig. 2-b. The distance between the 
indentation wall and the center of these vortices is measured from the instantaneous streamlines and plotted 
versus time in Fig. 2-a. The dashed lines represent the present computations. The solid lines denote 
computational results from the fractional step approach implemented by Rosenfcld.^ 3 Dotted lines are numerical 
results of stream function vorticity formulations from Ralph and Pedley. 18 Experimental measurements by Pedley 
and Stephanoff 17 are represented by square symbols. The agreement between numerical and experimental 
results is fairly good. There is a discrepancy between numerical results and experimental measurements about 
the location of the vortex A. The present results and Rosenfeld’s computations predict the location of vortex A 
0.4 units closer to the indentation than experimental findings. However, the locations of the remaining vortices 
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are correctly predicted if the distance is measured from the center of the vortex A. This underprediction of the 
separation length of vortex A is thought to be caused by an inaccurate description of the indentation wall shape. 
The present study and Rosenfeld’s study used the same grid and the same wall shape. Even though the solution 
algorithms of the two computational studies are completely different, the agreement between the numerical 
results is good. For these reasons, only the location of vortices B, C, and D are compared with the experimental 
measurements. 


Flow Through Tilting Disk Valve 

This problem was chosen to develop and validate a procedure which will be used for the valve region of the Penn 
State artificial heart. In the Bjork-Shiley tilting disk heart valve, the tilting disk is placed in front of the sinus 
region of the human aorta. The aortic root has three sinuses about 120 degrees apart from one another. The 
tilting disk valve model used in this computation is simplified by assuming that the sinus region of the aorta has 
a circular cross-section. The cage and struts which hold the free-floating disk inside of the sewing ring are not 
included in the geometry. It is also assumed that the walls do not have an elastic deformation. The channel 
length is taken to be five aorta diameters long. The computational geometry used in these unsteady flow 
computations is given in Fig. 3. The disk motion is illustrated by showing three different positions of the disk, 
at angles of 75, 50, and 30 degrees as measured from the centerline of the aorta. The tilting disk is allowed to 
rotate about the horizontal axis that is 1/6 of a disk diameter below the center of the disk. Because of this 
asymmetric disk orientation, the flow is three dimensional. 

The Chimera grid embedding technique, which has been successfully used for external flow problems, has been 
employed by using two overlapped grids as shown in Fig. 4. Grid 1 occupies the whole region in the aorta from 
entrance to exit, and remains stationary. Grid 2 wraps around the tilting disk, and moves with the disk. In the 
Chimera grid embedding technique, grid points which lie within the disk geometry and outside the channel grid 
are excluded from the solution process. These excluded points are called hole points, and the immediate 
neighbors of the hole points are called fringe points. The information is passed from one grid to another one 
via fringe and grid boundary points by interpolating the dependent variables. Tri-linear interpolation is used in 
the present computations. In order to distinguish the hole and fringe points from regular computational points, 
an IBLANK array is used in the flow solver. For hole, grid boundary, and fringe points IBLANK is set to zero, 
otherwise it is set to one. In order to exclude the hole and grid boundary points from the solution procedure, 
the coefficients of the system of algebraic equations and right-hand-side terms are multiplied by the IBLANK 
value. If the grid point is a hole, an outer boundary, or a fringe point the value of (1— IBLANK) is added to the 
main diagonal of the matrix equation. 

Presented here are the results of steady flow with a fixed disk angle and unsteady flow with the disk motion in 
the configuration described above. The problems are nondimensionalized by using the entrance diameter as a 
unit length, and the average inflow velocity as the unit velocity. In order to reduce the computational effort and 
memory size, the inflow and outflow boundaries are placed a short distance from the region of interest in 
comparison with the boundaries in the experimental studies. In addition, the exact shape of the sinus region of 
the aorta used in the experiments is not known. These discrepancies could lead to slight differences between 
present computations and experimental measurements. 

Steady-state calculations for the 30 degree disk orientation have been carried out for Reynolds numbers in the 
range of 2000 to 6000, in which experimental data are available. The Reynolds number is based on the diameter 
and the mean velocity at the entrance of the channel. A mixing length algebraic turbulence model which is used 
for incompressible flow through the space shuttle main engine turnaround duct is utilized. Figure 5 shows the 
pressure drop across the Bjork-Shiley tilting disk valve at different flow rates of physiological interest. The 
computed and measured axial velocity profiles at 42 mm downstream from the disk are shown in Fig. 6. Axial 
velocity profiles are plotted in the horizontal plane through the center of the channel. The numerical results are 
shown with dots and the experimental results are shown with triangles. The numerical results compare favorably 
with the experimental measurements. 3 The largest discrepancy is seen near the walls, where the boundary layer 
is overestimated by the calculation. Figure 7 shows the velocity vectors at five longitudinal stations. The flow, 
which is directed to the upper part of the aorta, generates vortices in the sinus region of the aorta and a large 
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separated region along the lower wall of the aorta. Since separated and low flow regions have potential for 
thrombus formation, clotting may occur on the upper sinus region and the lower wall of the aorta. Figliola and 
Mueller 5 also present mean velocity profiles, which show similar flow characteristics to those indicated by by 
computational study, at several locations. They computed the shear stress from the measured velocity field and 
observed that the maximum shear occurs at the top wall downstream of the sinus region of the aorta. This is in 
agreement with the velocity plot shown in Fig. 7, in which there are large velocity components just off the wall 
in that location. Particle traces in Fig. 8 indicate that the flow does not separate adjacent to the tilting disk;. The 
tilting disk separates the flow into a major flow region, which is along the upper wall of the tube, and a minor 
flow region along the lower wall of the tube. Separation, reverse flow and swirling motion mostly occur in the 
minor flow region. 

Figure 9 shows vorticity magnitude contours on the surface of the tube, outflow surface, and inflow surface of 
the disk, respectively. It is assumed that maximum vorticity magnitudes indicate the regions of high shear. The 
sewing ring surface and the edges of the disk are the regions having maximum vorticity magnitude. The upper 
wall of the channel also has considerably high vorticity magnitudes. 

Unsteady flow calculations have been carried out in order to demonstrate and analyze the flow during disk 
opening and closing. For the present computation, one cycle of valve opening and closing requires 70 physical 
time steps. During each time step, subiterations are carried out until the maximum divergence of velocity and 
maximum residual drop below 10." 3 The computing time required for one cycle of the valve opening and closing 
is approximately 5 Cray-YMP hours. During the valve opening, inflow velocity is imposed at the entrance of the 
channel. The inflow velocity is chosen as a sine function in time. The disk rotation is specified as a linear function 
in time. Since the forces acting on the disk are known from the numerical solution, the disk rotation angle can 
be determined. However, there is a limitation in the disk rotation angle. For large disk rotation angle, some 
information may be lost between the grids when the grid embedding technique is used. In order to prevent the 
information loss, the maximum allowed disk rotation angle at each physical time step is taken to be less than 
three degrees. 

Figures lO-a through lO-c illustrate the velocity vectors on the lateral symmetry plane at t/T = 0.128, 0.257, and 
0.385 respectively. The velocities are very high in the region between the disk and the channel wall as shown in 
Fig. lO-a. During the disk opening, two vortices are formed at the upper and lower edges of the disk. The flow 
starts to separate behind the disk and reattaches to the wall as shown in Fig. lO-b. The stagnation region behind 
the disk moves downstream as the disk rotates. Highly skewed velocity profiles are seen downstream from the 
disk as illustrated in Fig. lO-c. The growth of the vortices has also been observed in the sinus region of the aorta 
while the flow opens the valve. Along the lower wall a separation region is formed. 

Artificial Heart Flow 

The geometry of the Penn State artificial heart model is composed of a cylindrical chamber with two tube 
extensions (see Fig. 11). The inflow (mitral) and outflow (aortic) tubes contain concave tilting disks which open 
and close to act as valves. In the computational model, tilting disk mitral valve orientation in time was obtained 
from the experimental data provided by the Penn State University. The aortic valve orientation in time was 
approximated to mitral valve orientation with a phase difference. The pumping action is provided by a pusher 
plate whose velocity is sinusoidal in time. Pusher plate diameter is 7.26 cm, with a stroke length of 2.28 cm. The 
problem is nondimensionalized with the inflow tube diameter, which is 2.54 cm, and a unit velocity of 20 cm /sec. 
In the computational study, the Reynolds number based on the unit length and velocity is 900. Initially, the flow 
was started at rest, and four cycles of the pumping action were completed using a Cray-YMP computer at 
NASA-Ames Research Center. One cycle of the pusher plated motion required 240 physical time steps. At each 
time step, the equations were iterated until the maximum divergence of velocity was reduced below 10.' 2 During 
most of the cycle 10-20 subiterations were required (for more detail, see Refs. 14-15). 

In order to handle the geometric complexity and the moving boundary problems, a zonal method and an 
overlapped grid embedding 16 scheme are employed, respectively. In the zonal method, a complex computational 
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domain is divided into several simple subdomains. The overlapped grid embedding scheme allows subdomains 
to move relative to each other, and provides great flexibility when the boundary movement creates large 
displacements. The computational grid for this heart model is shown in Fig. 11. Grid 1 is generated for the 
pusher plate and moves with it. Grid 2 occupies the chamber and remains stationary. Grid 3 and Grid 5 are for 
the inflow and outflow tube extensions, respectively. Grid points for the tubes and grid points for the chamber 
are overlapped on three common planes. In other words, the grid points for the tubes start three stencils inside 
of the chamber outer boundaries. Zonal boundary conditions are used at the interface boundaries. Grid 4 and 
Grid 6 wrap around tilting disks, and move with the disks. An overlapped grid embedding scheme is employed 
between moving grids and stationary grids. 

Computed results presented next are mainly to demonstrate how CFD can be used to understand the flow in 
the artificial heart, and comparison with experiment is qualitative. Unsteady particle traces are illustrated in Fig. 

12. The particles are released near the inflow valve at the beginning of the fourth cycle. The figure is plotted at 
non-dimensional time t/T - 0.45 into the period at which time the pusher plate is close to its lowest position, 
where T denotes the period for the pusher plate’s motion. Figure 12 shows that the flow creates a strong vortex 
in the center region of the chamber. The particles have a swirling motion against the back wall opposite the 
mitral valve opening. The flow also separates at the connection region of the chamber and the inflow tube. Figure 
13 shows the computed velocity vectors on the horizontal mid-plane at non-dimensional time 0.375. At that time 
the pusher plate is moving down, and the mitral valve is opened. Figure 13 also indicates the presence of strong 
circulation in the chamber. However, the three dimensional structure of the flow can not be seen clearly because 
the vectors are plotted on a two dimensional plane. 

The strong vortex in the center of the chamber is actually created where the chamber and inflow tube are 
connected. The vortex moves to the core of the chamber in time. Experimental measurements by Baldwin and 
Tarbel 19 are illustrated in Fig. 14. Since the computational study does not include the blood sac inside the 
chamber, the comparison between experimental and computational results is qualitative. In addition, the 
Reynolds number in the experimental study is 1.7 times larger than the Reynolds number in the computational 
study because the flow is assumed laminar in the present computations. The biggest discrepancy between 
experimental measurements and computational results is the location of the vortex core in the chamber. In Fig. 

13, the vortex is off center in the chamber. In Fig. 14, the vortex is located almost in the center of the chamber. 
Another difference can be seen in the wake region of the mitral valve. Since the Reynolds number in the 
computational study is lower, the wake is not as strong as the wake in Fig. 14. The Reynolds number in the 
future computational study will be increased by including the turbulence modeling in order to have a quantitative 
comparison between experimental and computational results. 

During the second half of the cycle time, the pusher plate moves upward, and the outflow valve is opened. A top 
view of the computed velocity vectors on the horizontal plane at t/T = 0.625 is plotted in Fig. 15. Since the 
inflow valve is closed, residual eddies are quite large near the disk. However, they are quickly weakened as the 
pressure builds up inside the chamber. Measured velocity vectors at t/T = 0.625 are shown in Fig. 16. 

Figure 17 shows additional vortex structures in the vertical plane of the chamber at non-dimensional time 0.25. 
That vortex structure causes the swirling motion of the particles previously mentioned in Fig. 12. The swirling 
flow in the outflow tube is shown in Fig. 18. The velocity vectors in the cross-sectional plane of the outflow tube 
at t/T = 0.75 are plotted. The cross-sectional plane is one non-dimensional unit downstream of the tilting disk 
valve, and the normal vector of that cross-section is in positive x-direction shown in Fig. 13. 

SUMMARY 

An efficient and robust solution procedure is developed, and validated for numerical simulations of internal flows 
through artificial heart devices. The solution procedure for unsteady incompressible viscous flow computations 
has been extended with the incorporation of the grid embedding approach. This has been used to simulate the 
flow through a tilting disk heart valve and the flow through the Penn State artificial heart model. Separated and 
secondary flow regions have been pointed out in the tilting disk heart valve and artificial heart flow simulations. 
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The vortex created in the central portion of the Penn State artificial heart provides good wall washing over the 
entire chamber. The present capability of simulating complicated internal flow problems with moving boundaries 
is demonstrated. The procedure developed in this study is quite general and applicable for various types of 
artificial heart and valve geometries. It is hoped that researchers and designers in the artificial heart research 
may further benefit from the computational ability obtained in the current work. 
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